/*
 * 	This file is a modification of phylip's drawtree.c. It is also part 
 *  of MVD_GUI.
 *
 *  MVD_GUI is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MVD_GUI is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with MVD_GUI.  If not, see <http://www.gnu.org/licenses/>.
 *  This version copyright 2010 Desmond Schmidt
 */
/* In addition to the above, the following notice also applies */
/* Version 3.6.  Copyright (c) 1986-2004 by the University of Washington and
  Written by Joseph Felsenstein and Christopher A. Meacham.  Additional code
  written by Sean Lamont, Andrew Keefe, Hisashi Horino, Akiko Fuseki, Doug
  Buxton and Michal Palczewski.  Permission is granted to copy, distribute,
  and modify this program provided that (1) this copyright message is
  not removed and (2) no fee is charged for this program. */

#include <string.h>
#include <stdio.h>
#include <math.h>
#include <syslog.h>
#include <stdlib.h>
#include "phylip.h"
#include "draw.h"
#include "drawtree.h"

#define GAP 0.5
#define MAXITERATIONS 100
#define MINIMUMCHANGE   0.0001

#ifndef INFINITY
#define INFINITY 0x7fffffff
#endif

FILE *plotfile;
FILE *intree;
char plotFilePath[FILE_NAME_LEN];
char inTreePath[FILE_NAME_LEN];
int maxNumOfIter;
int numlines;
int nextnode;
node *root;
node *grbg;
node *where;
pointarray nodep;
boolean goteof;
boolean haslengths;
boolean rotate;
boolean improve;
boolean nbody;
boolean regular;
boolean rescaled;
boolean labelavoid;
boolean uselengths;
double topoflabels;
double bottomoflabels;
double rightoflabels;
double leftoflabels;
double maxx;
double maxy;
double minx;
double miny;
double xmargin;
double xsize;
double bscale;
double ymargin;
double ysize;
double maxchange;
double xcorner;
double ycorner;
double xunitspercm;
double yunitspercm;
double labelheight;
double charht;
double yscale;
double expand;
double pagex;
double pagey;
double paperx;
double papery;
double hpmargin;
double vpmargin;
double xscale;
double xoffset;
double yoffset;
double treeangle;
double ark;
double labelrotation;
double *textlength, *firstlet;
plottertype plotter;
labelorient labeldirec;
growth grows;
char fontname[FONT_NAME_LEN];
/* calculate force acting between 2 nodes and return the force in pForce. 
   Remember to pass the index subnodes to this function if needed. 
   Force should always be positive for repelling.  Angle changes to 
   indicate the direction of the force.  The value of INFINITY is the cap
   to the value of Force.  
   There might have problem (error msg.) if pFromSubNode and pToSubNode 
   are the same node or the coordinates are identical even with double 
   precision.  */
void force_1to1(node *pFromSubNode, node *pToSubNode, double *pForce,
	double *pAngle, double medianDistance)
{
  	double distanceX, distanceY, distance, norminalDistance;
  	distanceX = pFromSubNode->xcoord - pToSubNode->xcoord;
  	distanceY = pFromSubNode->ycoord - pToSubNode->ycoord;
  	distance = sqrt( distanceX*distanceX + distanceY*distanceY );
  	norminalDistance = distance/medianDistance;

  	if (norminalDistance < epsilon)
  	{
    	*pForce = INFINITY;
  	}
  	else
  	{
    	*pForce = (double)1 / (norminalDistance * norminalDistance);
    	if (*pForce > INFINITY) *pForce = INFINITY;
  	}
  	*pAngle = computeAngle(pFromSubNode->xcoord, pFromSubNode->ycoord,
    	pToSubNode->xcoord, pToSubNode->ycoord);
  return;
}
/* pToSubNode is where all the relevent nodes apply forces to.
   All branches are visited except the branch contains pToSubNode.
   pToSubNode must be one of the branch out of the current Node (= out of one
   of the subnode in the current subnodes set.)
   Most likely pPivotSubNode is not the index subNode!  In any case, only
   the leafs are consider in repelling force; so, no worry about index subNode.
   pTotalForce and pAngle must be set to 0 before calling this function for the
   first time, or the result will be invalid.  pPivotSubNode is named for 
   external interface.  When calling totalForceOnNode() recursively, 
   pPivotSubNode should be thought of as pFromSubNode.
 */
void totalForceOnNode(node *pPivotSubNode, node *pToSubNode,
	double *pTotalForce, double *pAngle, double medianDistance)
{
  	node *pSubNode;
  	double force, angle, forceX, forceY, prevForceX, prevForceY;

  	pSubNode = pPivotSubNode;

  	/* visit the rest of the branches of current node; the branch attaches to 
     	the current subNode may be visited in the code down below. */
  	while (pSubNode->next != NULL && pSubNode->next != pPivotSubNode)
  	{
    	pSubNode = pSubNode->next;
    	if ( pSubNode->back != NULL && pSubNode->back != pToSubNode) 
      		totalForceOnNode(pSubNode->back, pToSubNode, pTotalForce, pAngle, 
          	medianDistance);
  	}

  	/* visit this branch; You need to visit it for the first time - at root only!
   	 *
     * Modified so that all nodes are visited and calculated forces, instead of 
     * just the leafs only. 
     * use pPivotSubNode instead of pSubNode here because pSubNode stop short
     *  just before pPivotSubNode (the entry node) */
  	if ( pPivotSubNode == root && pPivotSubNode->back != NULL 
       	&& pPivotSubNode->back != pToSubNode) 
    totalForceOnNode(pPivotSubNode->back, pToSubNode, pTotalForce, pAngle, 
        medianDistance);
  	/* Break down the previous sum of forces to components form */
  	prevForceX = *pTotalForce * cos(*pAngle);
  	prevForceY = *pTotalForce * sin(*pAngle);
  	force_1to1(nodep[pPivotSubNode->index-1], pToSubNode, &force, &angle, 
      medianDistance);
  	/* force between 2 nodes */
  	forceX = force * cos(angle);
  	forceY = force * sin(angle);
  	/* Combined force */
  	forceX = forceX + prevForceX;
  	forceY = forceY + prevForceY;
  	/* Write to output parameters */
  	*pTotalForce = sqrt( forceX*forceX + forceY*forceY );
  	*pAngle = computeAngle((double)0, (double)0, forceX, forceY);
  	return;
}
/* Return the equivalent value of angle that is within
   0 to 2*pi */
double capedAngle(double angle)
{
  	while (angle < 0 || angle >= 2*pi)
  	{
    	if(angle < 0)
    	{
      		angle = angle + 2*pi;
    	}
    	else if (angle >= 2*pi)
    	{
      		angle = angle - 2*pi;
    	}
  	}
  	return angle;
}
/* it return the sign of the moment caused by the force, applied
   to the tip of the refereceVector; the root of the refereceVector
   is the pivot. */
double signOfMoment(double xReferenceVector, double yReferenceVector,
                double xForce, double yForce)
{
  	double angleReference, angleForce, sign;

  	angleReference = computeAngle((double)0, (double)0, xReferenceVector, 
                                yReferenceVector);
  	angleForce = computeAngle((double)0, (double)0, xForce, yForce);
  	angleForce = capedAngle(angleForce);
  	angleReference = capedAngle(angleReference);

  	/* reduce angleReference to 0 */
  	angleForce = angleForce - angleReference;
  	angleForce = capedAngle(angleForce);

  	if (angleForce > 0 && angleForce < pi)
  	{
    	/* positive sign - force pointing toward the left of the reference 
       	line/vector.  */
    	sign = 1;
  	}
  	else
  	{
    	/* negative sign */
    	sign = -1;
  	}
  	return sign;
}
/* added - danieyek 990205 */
/* pStackTop must be the current top element of the stack, where we pop an
   element from the top of it.  
   ppStackTop must be the location where we can find pStackTop.
   This function "returns" the revised top (element) of the stack through 
   the output parameter, ppStackTop. 
   The last element on the stack has the "back" (pStackElemBack) pointer
   set to NULL.  So, when the last element is poped, ppStackTop will be
   automatically set to NULL.  If popNodeFromStack() is called with 
   ppStackTop = NULL, we assume that it is the error caused by over popping
   the stack.
*/
void popNodeFromStack(stackElemType **ppStackTop, node **ppNode)
{
  	stackElemType *pStackT;

  	if (ppStackTop == NULL)
  	{
    	syslog(LOG_ERR,"drawtree - a call to pop while the stack is empty.\n");
    	exit(1);
  	}
  	pStackT = *ppStackTop;
  	*ppStackTop = pStackT->pStackElemBack;
  	*ppNode  = pStackT->pNode;
  	free(pStackT);
  	return;
}
void polartrav(node *p, double xx, double yy, double firstx,
    double firsty, double *leftx, double *lefty, double *rightx,
	double *righty)
{
  	/* go through subtree getting left and right vectors */
  	double x, y, xxx, yyy, labangle = 0;
  	boolean lookatit;
  	node *pp;

  	lookatit = true;
  	if (!p->tip)
    	lookatit = (p->next->next != p || p->index != root->index);
  	if (lookatit) 
	{
    	x = nodep[p->index - 1]->xcoord;
    	y = nodep[p->index - 1]->ycoord;
    	if (p->tip) 
		{
      		if (labeldirec == fixed) 
			{
        		labangle = pi * labelrotation / 180.0;
        		if (cos(p->oldtheta) < 0.0)
          			labangle = labangle - pi;
      		}
      		if (labeldirec == radial)
        		labangle = p->theta;
      		else if (labeldirec == along)
        		labangle = p->oldtheta;
      		else if (labeldirec == middle)
        		labangle = 0.0;
      		xxx = x;
      		yyy = y;
      		if (labelavoid) 
			{
        		if (labeldirec == middle) 
				{
          			xxx += GAP * labelheight * cos(p->oldtheta);
          			yyy += GAP * labelheight * sin(p->oldtheta);
          			xxx += labelheight * cos(labangle) * textlength[p->index - 1];
          			if (textlength[p->index - 1] * sin(p->oldtheta) < 1.0)
            			xxx += labelheight * cos(labangle) * textlength[p->index - 1];
          			else
            			xxx += 0.5 * labelheight * cos(labangle)
                     		* textlength[p->index - 1];
          			yyy += labelheight * sin(labangle) * textlength[p->index - 1];
        		}
        		else 
				{
          			xxx += GAP * labelheight * cos(p->oldtheta);
          			yyy += GAP * labelheight * sin(p->oldtheta);
          			xxx -= labelheight * cos(labangle) * 0.5 * firstlet[p->index - 1];
          			yyy -= labelheight * sin(labangle) * 0.5 * firstlet[p->index - 1];
          			xxx += labelheight * cos(labangle) * textlength[p->index - 1];
          			yyy += labelheight * sin(labangle) * textlength[p->index - 1];
        		}
      		}
      		if ((yyy - yy) * firstx - (xxx - xx) * firsty < 0.0) 
			{
        		if ((yyy - yy) * (*rightx) - (xxx - xx) * (*righty) < 0.0) 
				{
          			(*rightx) = xxx - xx;
          			(*righty) = yyy - yy;
        		}
      		}
      		if ((yyy - yy) * firstx - (xxx - xx) * firsty > 0.0) 
			{
        		if ((yyy - yy) * (*leftx) - (xxx - xx) * (*lefty) > 0.0) 
				{
          			(*leftx) = xxx - xx;
          			(*lefty) = yyy - yy;
        		}
      		}
    	}
    	if ((y - yy) * firstx - (x - xx) * firsty < 0.0) 
		{
      		if ((y - yy) * (*rightx) - (x - xx) * (*righty) < 0.0) 
			{
        		(*rightx) = x - xx;
        		(*righty) = y - yy;
      		}
    	}
    	if ((y - yy) * firstx - (x - xx) * firsty > 0.0) 
		{
      		if ((y - yy) * (*leftx) - (x - xx) * (*lefty) > 0.0) 
			{
        		(*leftx) = x - xx;
        		(*lefty) = y - yy;
      		}
    	}
  	}
  	if (p->tip)
    	return;
  	pp = p->next;
  	while (pp != p) 
	{
    	if (pp != NULL)
      		polartrav(pp->back,xx,yy,firstx,firsty,leftx,lefty,rightx,righty);
    	pp = pp->next;
  	}
}
void polarize(node *p, double *xx, double *yy)
{
  	double TEMP, TEMP1;

  	if (fabs(p->xcoord - (*xx)) > epsilon)
    	p->oldtheta = atan((p->ycoord - (*yy)) / (p->xcoord - (*xx)));
  	else if (p->ycoord - (*yy) > epsilon)
      	p->oldtheta = pi / 2;
  	if (p->xcoord - (*xx) < -epsilon)
    	p->oldtheta += pi;
  	if (fabs(p->xcoord - root->xcoord) > epsilon)
    	p->theta = atan((p->ycoord - root->ycoord) / (p->xcoord - root->xcoord));
  	else if (p->ycoord - root->ycoord > 0.0)
      	p->theta = pi / 2;
    else
      	p->theta = 1.5 * pi;
  	if (p->xcoord - root->xcoord < -epsilon)
    	p->theta += pi;
  	TEMP = p->xcoord - root->xcoord;
  	TEMP1 = p->ycoord - root->ycoord;
  	p->r = sqrt(TEMP * TEMP + TEMP1 * TEMP1);
}
/* added - danieyek 990128 */
/* After calling tilttrav(), if you don't polarize all the nodes on the branch
   to convert the x-y coordinates to theta and radius, you won't get result on 
   the plot!  This function takes a subnode and branch out of all other subnode
   except the starting subnode (where the parent is), thus converting the x-y 
   to polar coordinates for the branch only.  xx and yy are purely "inherited"
   features of polarize().  They should have been passed as values not 
   addresses. */
void polarizeABranch(node *pStartingSubNode, double *xx, double *yy)
{
  	node *pSubNode;

  	pSubNode = pStartingSubNode;
  	/* convert the current node (note: not subnode) to polar coordinates. */
  	polarize( nodep[pStartingSubNode->index - 1], xx, yy);

  	/* visit the rest of the branches of current node */
  	while (pSubNode->next != NULL && pSubNode->next != pStartingSubNode)
  	{
    	pSubNode = pSubNode->next;
    	if ( pSubNode->tip != true )
    	  	polarizeABranch(pSubNode->back, xx, yy);
  	}
  	return;
}
/* traverse to move successive nodes */
void tilttrav(node *q, double *xx, double *yy, double *sinphi,
	double *cosphi)
{
  	double x, y;
  	node *pp;

  	pp = nodep[q->index - 1];
  	x = pp->xcoord;
  	y = pp->ycoord;
  	pp->xcoord = (*xx) + (x - (*xx)) * (*cosphi) + ((*yy) - y) * (*sinphi);
  	pp->ycoord = (*yy) + (x - (*xx)) * (*sinphi) + (y - (*yy)) * (*cosphi);
  	if (q->tip)
    	return;
  	pp = q->next;
  	while (pp != q) 
	{
	    /*
	    if (pp != root)
	    */
	    if (pp->back != NULL)
 	     	tilttrav(pp->back,xx,yy,sinphi,cosphi);
    	pp = pp->next;
  	}
}
/* added - danieyek 990204 */
/* pStackTop must be the current top element of the stack, where we add another
   element on top of it.  
   ppStackTop must be the location where we can find pStackTop.
   This function "returns" the revised top (element) of the stack through 
   the output parameter, ppStackTop. 
   The last element on the stack has the "back" (pStackElemBack) pointer
   set to NULL.  So, when the last element is poped, ppStackTop will be
   automatically set to NULL.  If popNodeFromStack() is called with 
   ppStackTop = NULL, we assume that it is the error caused by over popping
   the stack.
*/
void pushNodeToStack(stackElemType **ppStackTop, node *pNode)
{
  	stackElemType *pStackElem;

  	if (ppStackTop == NULL)
  	{
    	/* NULL can be stored in the location, but the location itself can't 
       	be NULL! */
    	syslog(LOG_ERR,"drawtree - error using pushNodeToStack(); " 
        "ppStackTop is NULL.\n");
    	exit(-1);
  	}
  	pStackElem = (stackElemType*)Malloc( sizeof(stackElemType) );
  	pStackElem->pStackElemBack = *ppStackTop;
  	/* push an element onto the stack */
  	pStackElem->pNode = pNode;
  	*ppStackTop = pStackElem;
  	return;
}
/* polar coordinates of a node relative to start */
void plrtrans(node *p, double theta, double lower, double upper)
{
  	int num;
  	double nn, pr, ptheta, angle, angle2, subangle, len;
  	node *pp, *qq;

  	nn = p->width;
  	subangle = (upper - lower) / nn;
  	qq = p;
  	pp = p->next;
  	if (p->tip)
    	return;
  	angle = upper;
  	do 
	{
	    angle -= pp->back->width / 2.0 * subangle;
	    pr = p->r;
	    ptheta = p->theta;
	    if (regular) 
		{
      		num = 1;
      		while (num * subangle < 2 * pi)
        	num *= 2;
      		if (angle >= 0.0)
        		angle2 = 2 * pi / num * (long)(num * angle / (2 * pi) + 0.5);
      		else
        		angle2 = 2 * pi / num * (long)(num * angle / (2 * pi) - 0.5);
   		} 
		else
      		angle2 = angle;
    	if (uselengths)
      		len = fabs(pp->back->oldlen);
    	else
      		len = 1.0;
    	pp->back->r = sqrt(len * len + pr * pr + 2 * len * pr * cos(angle2 - ptheta));
    	if (fabs(pr * cos(ptheta) + len * cos(angle2)) > epsilon)
      		pp->back->theta = atan((pr * sin(ptheta) + len * sin(angle2)) /
                 (pr * cos(ptheta) + len * cos(angle2)));
    	else if (pr * sin(ptheta) + len * sin(angle2) >= 0.0)
      		pp->back->theta = pi / 2;
    	else
      		pp->back->theta = 1.5 * pi;
    	if (pr * cos(ptheta) + len * cos(angle2) < -epsilon)
      		pp->back->theta += pi;
    	if (!pp->back->tip)
      		plrtrans(pp->back, pp->back->theta,
            	angle - pp->back->width * subangle / 2.0,
                angle + pp->back->width * subangle / 2.0);
    	else
      		pp->back->oldtheta = angle2;
    	angle -= pp->back->width / 2.0 * subangle;
    	pp = pp->next;
  	} while (((p != root) && (pp != qq)) || ((p == root) && (pp != p->next)));
}
